A Superglass Phase of Interacting Bosons 
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We introduce a Bose-Hubbard Hamiltonian with random disordered interactions as a model to 
study the interplay of superfluidity and glassiness in a system of three-dimensional hard-core bosons 
at half-filling. Solving the model using large-scale quantum Monte Carlo simulations, we show 
that these disordered interactions promote a stable superglass phase, where superflow and glassy 
density localization coexist in equilibrium without exhibiting phase separation. The robustness 
of the superglass phase is underlined by its existence in a replica mean-field calculation on the 
infinite-dimensional Hamiltonian. 



Despite the simplicity of its constituent atoms, the 
phase diagram of bulk helium has proven to be com- 
pelling and controversial - especially in light of recent 
observations of supersolid behavior in the crystal phase 
ofHelium-4 (^He) Several experiments have shown 

that superflow is enhanced by imperfections in the crystal 
lattice, such as roaming defects and interstitial atoms [3|. 
Most strikingly, new experiments by Davis and collabo- 
rators indicate that the onset of superflow is intimately 
tied to relaxation dynamics characteristic of glassy, or 
amorphous, solids [4]. The precise relationship between 
superflow and glassiness, and the possible existence of 
a novel superglass (SuG) state, is not without contro- 
versy in the experimental community (Hf. Even outside 
of the ^He context, the possible existence of a bosonic 
superglass phase is of broad theoretical importance - 
one that can be explored for example in simple models 
of interacting bosons. Lattice models have been instru- 
mental in explaining some of the basic phenomenology 
of strongly-coupled quantum systems, including that of 
^He. The prototypical model is the lattice Bose-Hubbard 
(BH) model ^ . Its extended phase diagram is known to 
contain superfluids, Mott-insulating "crystals", and su- 
persolid phases, where the latter refers to a state with co- 
existing superfluidity and broken translational symmetry 
in the particle density f?\. 

Models of interacting bosons with disorder have been 
considered for some time [6, 8] typically in the local chem- 
ical potential, as might be realized in ^He absorbed in 
porous media. In these and related models [9] arises 
either a "Bose-glass" (BG) phase, with localized disor- 
der in the particle densities but no coexisting superflow, 
or superfluid phases with locally inhomogeneous super- 
flow [lO[ - neither of which correspond to a SuG state. 
This can be understood in part by considering the na- 
ture of the states in the BG close to the Mott lobes 
which, through Anderson localization [ll|, essentially 



form single-particle localized states. For this reason, the 
BG cannot support phase coherence; other types of in- 
teractions beyond the random local chemical potential 
are necessary to induce superglassiness. In this paper we 
show that a thermodynamic SuG state of bosons can be 



stabilized via random pairwise boson-boson interactions. 

It is widely believed that disorder (in the form of dis- 
locations, grain boundaries, impurities, etc.) plays a role 
in the supersohd behavior of ^He although the pre- 
cise nature of this role is controversial. At the micro- 
scopic level, ^He atoms interact via a deterministic pair- 
wise potential in the continuum. If perfect crystal order 
is avoided by an effective quenching [l2[ into a glassy 
phase at low temperature, even without explicit random- 
ness in the Hamiltonian, both disorder and random frus- 
tration can be self-generated from the deterministic po- 
tential, similar to scenarios put forward in discussions of 
the structural glass transition [l3|. Therefore, we pro- 
pose that an effective theory to describe glassy behavior 
in ^He via this disordered and frustrated environment is 
a BH lattice model with randomly frustrated boson-boson 
interactions. We note that the phenomenology of random 
frustration and glassiness can also be constructed via a 
mathematical embodiment of a BH model on a random 
graph Q - an approach that leaves open the question of 
whether a SuG exists for real finite dimensional lattices. 

While the relevance of our model to the physics of ^He 
ultimately relies on microscopic details via which disor- 
der and frustration may be self-generated, it is of direct 
pertinence to collective phenomena in ultracold atomic 
gases. There, disorder may be introduced through mech- 
anisms such as speckle potentials or superimposed optical 
lattices with incommensurate periods [lij . 

Accordingly, we introduce the hard-core BH model, 

H = -Y.Vi,{n, - 1/2)K- - 1/2) -tY, {b\hj + 

(1) 

where {ij) indicates nearest-neighbor lattice sites. Here, 
rii is the density operator of bosons hopping on a three- 
dimensional (3D) cubic lattice, hi {h\) is the annihi- 
lation (creation) operator. Vij is a quenched random 
potential with bimodal distribution given by PiVij) = 
[5 {Vij —V)^5{yij + V)]/2, reminiscent of the canonical 
spin glass (SG) models studied extensively [l6|, [l^. The 
hopping (t) term is the standard kinetic energy term in 
the BH model, acting as a source of quantum fluctua- 



2 



tions which favor a superfiuid phase whenever V/t is suf- 
ficiently small. Eq. ([!]) is explicitly defined to maintain 
half-filling. Below, we use large-scale Stochastic Series 
Expansion (SSE) quantum Monte Carlo (QMC) simu- 
lations to demonstrate that these simple random disor- 
dered interactions can promote a stable superglass phase 
in the Hamiltonian ([T]), where superflow and glassy den- 
sity localization coexist in equilibrium without exhibiting 
phase separation. Following that, the existence of this 
superglass phase is corroborated by a replica mean-field 
calculation on an infinite-dimensional model. 
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FIG. 1: a) The superfiuid density, EA order parameter, and 
spin glass susceptibility as a function of V/t, for T/t = 1/4 
and L — 10. For V/t > 4.2, xsg suffers from algorithmic 
freezing or loss of ergodicity, making results unreliable with- 
out a refinement of simulation technique (e.g. parallel tem- 
pering), b) Finite-size scaling of the QMC observables in the 
SuG state at V/t = 4 and T/t = 1/4. In (a) and (b), statisti- 
cal error bars related to QMC sampling are smaller than the 
symbol size, c) Histogram of the three observables at L = 10 
in (b) using 500 RGBs. The x-axis has been divided by the 
mean of the distribution in each case, and the y-axis of each 
has been offset for clarity. There is no clear separation of the 
SuG RGBs into distinct glassy and superfiuid bins. 



QMC simulations are performed using the finite- 
temperature SSE technique [l8], modified to allow for 
the random disorder in Eq. ([1]) by employing different 
directed loop equations for each individual bond inter- 
action. A single simulation consists of constructing a 
3D cubic lattice (size N = L^) with some unique choice 
of specific bimodal bond distribution, Vij^ referred to as 
one realization of disorder (ROD). Each specific simula- 
tion, for a given system size, temperature, and interac- 
tion strength V/t, is equilibrated for 10^ MC steps before 
thermodynamic data is collected (for 2 x 10^ steps), this 
MC average denoted below by (...). The simulation pro- 
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FIG. 2: The phase diagram of the (a) quantum (Eq. [Tj) and 
(b) mean-field (Eq. [7|) models, solved with QMC and replica 
mean-field methods respectively. Both phase diagrams show a 
robust superglass (SuG) phase, as well as superfiuid (SF) and 
quantum spin glass (SG) states. The data points in a, used to 
determine the phase boundaries, were obtained when ps and 
^EA surpassed an approximate numerical value of 1/ on a 
A/'=12xl2xl2 system. It is interesting to note that the 
main qualitative difference between the two phase diagrams 
is the absence of a direct superfiuid to superglass transition 
in the mean- field solution, (b), as a function of temperature. 



cedure is repeated for many RODs - typically between 
10^ and 10^, and denoted as [...]avg- 

We collect data for thermodynamic quantities which 
are able to quantify both the superfiuid and the glassy 
characteristics of various phases of the model. Our defi- 
nition for superflow is based on the superfiuid density. 
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where Wk is the winding number [l9| measured in the 
lattice direction k. Two quantities commonly used in 
SG models to indicate a glassy phase are the Edwards- 
Anderson (EA) order parameter. 
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and the replica overlap parameter [l6| . 

qah = ^^{n.-l/2r{n.-l/2)'', (4) 



where two identical (independently simulated) replicas 
of the system are labeled a and b for a single ROD. The 
disorder- averaged SG susceptibility is thus defined as 



XSG = N [{ql)] 



avg 



(5) 



Figure [T] a shows these observables as a function of V/t. 
For some finite region of V/t, the model displays robust 
superflow (characterized by finite ps) and a finite EA or- 
der parameter and SG susceptibility - these trends sur- 
vive into the thermodynamic limit (see Fig. [T]b). This 



indicates that the SuG is a stable equihbrium phase in 
this model, persisting down to zero temperature (Fig.[2j). 

The simplicity of Hamiltonian ([1]) suggests a clear 
mechanism for the formation of the SuG phase; it arises 
from the competition of disorder in the boson interactions 
with quantum fluctuations. One may therefore wonder 
whether this mechanism is robust enough to be found in 
standard analytical treatments of well-studied spin glass 
models. We address this question by employing an ex- 
act mapping of the hard-core boson model in Eq. ([T]) 
to a quantum spin model through the transformation 
to spin- 1/2 operators: Sf = 71^ — 1/2, S~ = 6^, and 
= b]. Equation [T] can then be written as a standard 
XXZ model: 



H = 



E 



(6) 



(ij) 



where Jxy = 2t and the off-diagonal part of the Hamilto- 
nian is the XY model, which (e.g. for V = 0) has in-plane 
ferromagnetic long-range order, {S^) 7^ 0, corresponding 
to the superfluid state in the boson language. One now 
recognizes the diagonal part of (j6j) as the standard EA 
Ising spin model, for which the expected SG phase is 
reproduced for our model as t ^ (Figs [1] and [2]). 

An analytical solution for the 3D EA model {t = 
0) is not available. However, the mean-field (MP) 
Sherrington-Kirkpatrick model can be solved exactl y by 
Parisi's replica symmetry breaking (RSB) ansatz [20|. 
Here, we generalize the replica MP method to include 
the effects from both the conventional Vij Ising coupling, 
and the off-diagonal Jxy-coup\mg [21]. The MP model is 

HMF = -Y.V,jS!S^-^Y.(^^^! + ^f^')^ (7) 

where the summation is over all distinct pairs of spins. 
Por simplicity, we use a Gaussian distribution, PiVij) = 

^ "^"^ ( 2V^'^ ) "> which does not qualitatively 
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change the physics compared to a bimodal PiVij). The 
disorder averaged free energy using the replica formalism 
is given by pN f = -[lnZ]avg = - limn^o ^([^""lavg - 1), 



where 



Javg 



denotes disorder averaging. We introduce 



the Hubbard- Stratonovich (HS) fields for SG overlap 
Q""^ = {S'^SD \f a^h, overlap within the same replica 
= {S^S^)^ and off-diagonal ordering parameter = 
{S^). We assume that W = M and W = R are replica 
independent and, for simplicity, consider solely a 1-step 
RSB for Q^^ . Following Parisi's parameterization for 1- 
step RSB [1^, we divide the n replicas into n/m groups 
of m replicas, and set Qah = Qi if a and h belong to the 
same group; or Qah = Qq \i a and h belong to different 
groups. Assuming that the HS fields are static, we obtain 
the free energy at 1-step RSB, 
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-m{Ql - Ql) 
2 j Dy cos\r{p/2^ HI + {JxyMf) 



,(8) 



where / Dx F{x) i 
and Hz = V{x^/Qq 



(i/v^)r 

-yVQi 



dx exp(-xV2)F(x) 



{pVf {R-Qi){l/2-R-Qi) 



Qo). The replica symmet- 
ric free energy can be recovered by setting any one of the 
conditions Qi = Qo, m = or m = 1. The self-consistent 
MP equations are obtained by optimizing the free energy. 
We find four phases analogous to those in the hard-core 
boson model (see Pig. [2]). In particular, a SuG phase, de- 
fined by coexisting in-plane order M 7^ (the magnetic 
analog to the superfluid order parameter), and SG over- 
lap parameters Qo ^ and Qi 7^ is clearly obtained. 
The vertical nature of the SuG-SF phase boundary orig- 
inates from the observation that the replica theory does 
not admit a self-consistent set of equations with nonzero 
M and Qo when V/t < 1. Physically, this results from 
the absence of phase fluctuations in the superfluid phase 
when the SG order parameter Qo is nonzero. 

In conclusion, both large-scale QMC simulations and 
replica MP calculations exhibit a robust, stable, equilib- 
rium superglass phase in the BH model, Eq. ([T]). The 
simplicity of the model reveals that the SuG phase de- 
velops under competition between random disordered 
interactions and quantum fluctuations. QMC simula- 
tions allow us to characterize the SuG phase in great 
detail; for example, as evident from the single ROD il- 
lustrated in Pig. m there is no phase separation into dis- 
tinct "super" and "glassy" regions within the simulation 
cell. Rather, the bulk of the particles in the system ex- 
hibit both off-diagonal superfluid and diagonal glassy or- 
dering. Most importantly, our results clearly indicate 
that glassiness can co-exist with superfluidity in a gen- 
eral bosonic system in three dimensions. Although we 
did not directly study a microscopic model of ^He, in the 
event that a metastable SuG [l2| does exist in ^He within 
experimentally-accessible timescales Q , quantum fluctu- 
ations and random frustration are the most relevant, if 
not inevitable, physical ingredients at play. Our model 
provides an affirmative answer that a SuG is possible un- 
der these two conditions. Furthermore, since the feasibil- 
ity of obtaining disordered interactions similar to our V 
term [Eq. ([T])] has recently been demonstrated in Fermi- 
Bose mixtures in random optical lattices 0, such cold 
atomic systems might also be suitable places to look for 
superglass phases in the future. With recent advances in 
both cold atom experiments and numerical methods [22| , 
the broad class of disordered Hamiltonian introduced in 
this paper provides new experimental and theoretical mo- 
tivation to study the interplay between superfluidity and 
quantum glassiness in numerous extensions of our model. 
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FIG. 3: The QMC simulation cell, with one random realization of disorder, averaged over 2 x 10^ Monte Carlo (MC) steps. 
Grey bonds correspond to Vij = V, striped bonds to Vij = —V. On the sites i of the lattice, the density localization is illustrated 
by the squared average of Sf = rii — 1/2, which is black for a site where a boson is delocalized (a "superfluid site"), white for 
a localized (frozen) site, and colored for partially localized Boson densities (0 < (iSf )^ < 0.25). The bond-strength, defined as 
the boson hopping operator (Bij) = {h\hj + hih^-) is represented as the thickness of the bonds. At bottom, the histograms of 
{Stf (left), and {Bij) (right) for each case: superfluid (a) dXV/t — 3, superglass (b) at V/t — 4, and spin glass (c) at V/t — 5. 
The double-peaked {Bij) histogram in (a) simply results from the bimodal distribution of Vij. 
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